function state = rk4(statep, e, h)
    k1 = h * dery(statep, e);
    k2 = h * dery(statep + 0.5 * k1, e);
    k3 = h * dery(statep + 0.5 * k2, e);
    k4 = h * dery(statep + k3, e);
    state = statep + (k1 + 2 * k2 + 2 * k3 + k4) / 6;
end

function dstate = dery(state, e)
    mass = 104305;
    S = 391.22;
    R0 = 6378137;
    g0 = 9.81;
    Vc = sqrt(R0*g0);
    rho0 = 1.2256;
    Hs = 7254.24;
    rad = state(1);
    lat = state(3);
    fpa = state(4);
    azi = state(5);
    V = sqrt(2*(1./rad - e));
    rho = rho0*exp(-(rad-1)*R0/Hs);
    if V*Vc > 4570
        alphadeg = 40;
    else
        alphadeg = 40 - 0.20705*(V*Vc - 4570)^2/340^2;
    end
    CL0 = -0.041065;
    CL1 = 0.016292;
    CL2 = 0.0002602;
    CD0 = 0.060505;
    CD1 = -0.03026;
    CD2 = 0.86495;
    CL = CL0 + CL1*alphadeg + CL2*alphadeg.^2;
    CD = CD0 + CD1*CL + CD2*CL.^2;
    L = 0.5*R0*rho.*V.^2*S.*CL/mass;
    D = 0.5*R0*rho.*V.^2*S.*CD/mass;
    raddot = sin(fpa)./D;
    londot = cos(fpa).*sin(azi)./(rad.*D.*cos(lat));
    latdot = cos(fpa).*cos(azi)./(rad.*D);
    fpadot = (L + (V.^2 - 1./rad).*cos(fpa)./rad)./(V.^2.*D);
    azidot = (V.^2.*cos(fpa).*sin(azi).*tan(lat)./rad)./(V.^2.*D);
    dstate = [raddot;londot;latdot;fpadot;azidot];
end